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Abstract. First-principle path integral Monte Carlo simulations were performed in order to 
analyze correlation effects in complex electron-hole plasmas, particularly with regard to the 
appearance of excitonic bound states. Results are discussed in relation to exciton formation in 
unconventional semiconductors with large electron hole mass asymmetry. 



1. Introduction 

Strongly coupled Coulomb systems have been in the focus of recent investigations in many fields, 
including plasma, nuclear, condensed matter and astrophysics 0|21II2II]- In these systems the 
Coulomb energy is of the order of, or even exceeds the mean kinetic energy. Besides, quantum 
effects are of vital importance. The formation of Coulomb liquids and solids is only one of the 
multi-faceted phenomena that might show up. Just now the existence of Coulomb crystals has 
been studied for neutral plasmas containing (at least) two oppositely charged components 0. 
Most notably crystallization of holes in semiconductors was predicted to occur in materials with 
a sufficiently large hole to electron effective mass asymmetry of about 80. 

Since bound state formation turned out to be the relevant limiting factor for the appearance 
of such an exceptional quantum hole (Coulomb) crystal at low temperatures [5], the aim of 
the present contribution is to analyze the formation of excitons in mass asymmetric electron- 
hole plasmas from first principles. A theoretical approach that is well suited to address this 
problem is the direct path integral quantum Monte Carlo (DPIMC) method for the N-particle 
density operator [5]. The DPIMC technique avoids additional approximations, such as the fixed 
node and restricted path integral method. It has been successfully applied to treat both strong 
interaction and quantum effects [HI E] • 

The organization of the paper is as follows. In Sec. 2 we briefly describe the essence of the 
DPIMC method we employed. Section 3 presents our simulation results. There we will give 
a detailed discussion of typical many-particle configurations at low and high temperatures, as 
well as of electron-electron, hole-hole and electron-hole pair distribution functions and structure 
factors. We conclude in Sec. 3, relating our results to the heavily debated exciton formation in 
intermediate valent Tm[Se,Te] 



2. Direct path integral Monte Carlo approach 

The starting point for the description of a two-component (equilibrium) quantum plasma is the 
partition function given at inverse temperature (3 = 1/ksT by 

Z(N e ,N h ,V,P) = / dqp(q,a;(3), (1) 

a V 

where q = {q e ,Qh} and a = {a e ,ah} denote the spatial coordinates and spin degrees of freedom 
of iV e electrons and N h holes (N e = N h ), respectively, i.e. q a = {qi, a ■ ■ ■ %a ■ ■ ■ <lN a ,a} and 
a a = {o"i, a • • • <y\,a ■ ■ ■ &N a ,a} with a = e,h. Then the pair distribution functions and the charge 
structure factors can be expressed as ^U] 

9ab(r) = z ^' N ^ h J2j d( l 6 ( r ha ~ ?l,a) K r 2,b ~ Q2,b) p(q, o\ ff) (2) 

and 

Q /.n _ Io° drr2 [9ab(r) ~ 1] sin(fcr)/(Ax) 
ab[) ~ \^drr^ 9ab {r)-l)\ ' W 

where r = ri jQ — r^fi- 

Of course, in general the exact density matrix p(q, a; (3) of an interacting quantum many- 
particle systems is not known but can be constructed using a path integral approach: 
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Jj2dqMp( q W,a;P)= Jdq^...dq^ 
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Here = p (g^, gW; A/?) = (g^ 1 ) | e - A/3 ^|gW) and A/3 = /3/(n + 1). The Hamiltonian for 
the electron- hole plasma (H = K + U) contains kinetic energy (K) and Coulomb interaction 
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Figure 1. (Color online) Bead 
representation of electrons and 
holes. Here Ag = 2-kTi 2 (3 /m e , 
\\ e = 2TTh 2 A(3/m e , gg = 

Qi°e + ^A,e^ e , and a = a'. 
The holes have a similar beads 
representation, however Xh is 
yjm e /mh times smaller. 



energy (U = U ee + U^h + U e h) parts. In Eq. (JIJ the index i = 1 . . . n + 1 labels the high- 
temperature [(n + l)ksT] density matrices pW. Accordingly each particle is represented by 
(re + 1) beads, i.e. the whole configuration of the particles is represented by a 3(N e + Nh)(n + 1)~ 
dimensional vector q = {qfj, g^ e +1) , qfj, . . . , . . . q^ ] ; qf^ . . . q^h}- Further details 

of the high-temperature path integral representation are given in Ref |S] • 

Figure. ^ illustrates the representation of one (light) electron and one (heavy) hole. The 
circle around the electron beads symbolizes the region that mainly contributes to the partition 
function path integral. The size of this region is of the order of the thermal electron wavelength 
A e (T), while typical distances between electron beads are of the order of the electron wavelength 
taken at (n + l)-times higher temperatures. The same representation is valid for holes but due 
to the larger hole mass the characteristic length scales are smaller by an factor of \Jm e jmy l . 
The spin gives rise to the spin part of the density matrix (S) with exchange effects accounted 
for by the permutation operators P e and i\ acting on the electron and hole coordinates g( ra+1 ) 
and spin projections a'. The sum is over all permutations with parity Kp e and Kp h . 

In the evaluation of the expressions obtained so far, we make use of a Monte Carlo 
scheme ^2] with different types of Monte Carlo steps: Either electron or hole coordinates (qi, a ) or 

dimensionless individual electron or hole beads (6 „) were moved until convergence is reached. 
Periodic boundary conditions are applied to the basic Monte Carlo cell. The procedure has 
been extensively tested. For example, the comparison with the known analytical expressions for 
pressure and energy of an ideal Fermi gas showed that Fermi statistics is well reproduced jSJ. 
Moreover we applied the method to few-electron systems in harmonic traps where again the 
analytically known limiting behavior of the energy was recovered [S]. For the present simulations 
of an mass-asymmetric electron-hole plasma, we have varied both particle and beads numbers 
and found that in order to obtain convergent results particle numbers N e = = 50 and 
beads numbers in the range of iVj = 20 are sufficient. Thereby the density matrix in the high 
temperature decomposition always was taken at temperatures above the electron-hole binding 
energy. Finally let us emphasize that the maximum statistical error of the results presented 
below is about 5% and can be systematically reduced by increasing the length of the Monte 
Carlo run. 

3. Simulation results 

We now use the DPIMC scheme outlined in the preceding section in order to study exciton 
formation in (non-ideal) electron- hole plasmas. Independent model parameters are the 
temperature (T), the masses of electrons (m e ) and holes (m/i), and the static dielectric constant 
(e). The plasma density is characterized by the Brueckner parameter r s , defined as the ratio of 
the mean distance between the particles d = [3/(47r(n e + n^))] 1 / 3 and the Bohr radius ag [TT] . 
where n e and are electron and hole densities, respectively. To make contact with experiments 
on the Tm[Te,Se] system we choose m e = 2, rrih = 80, and e = 25 (the masses are in units of 
the free electron mass). Then, assuming a simple Wannier exciton picture, the exciton binding 
energy would be about 500 K. 

3.1. Particle configurations 

Figure El shows typical "snapshots" of an electron-hole many-particle state with r s = 10 at low 
(50iT) and high (200-ftT) temperatures. Here, as a result of the temperature decomposition of 
the density matrix, each electron and hole is represented by several beads. Since the number of 
beads is twenty the high temperature density matrices correspond to a temperature of at least 
1000-fT, being about two times larger than the characteristic (excitonic) binding energy {E b x ). 

The spatial distribution of the beads reflects the position probability amplitude of the many- 
particle wave function, where the characteristic size of a certain cloud of beads is of the order 



Figure 2. (Color online) Snapshots of the electron-hole plasma at T = 50K (left) and T = 200K 
(right) (for r s = 10). Clouds of blue (yellow) dots represent electrons with spin up (down), 
whereas red (pink) dots mark the holes. Grey lines indicate the main simulation box which is 
periodically repeated in all spatial directions. 



of the thermal wave length of a single quantum particle. Because of the mass asymmetry the 
typical size of the electron clouds is about 6 times larger than that of the holes. Also, increase of 
the temperature leads to a reduce of the particle extension. Interestingly, at low temperatures 
practically all holes are closely covered by electron beads. This means electrons and holes form 
predominantly bound states and the whole system consists mainly of excitons. From Fig. |2] we 
can see that the average distance between particles (unbound electrons, holes and excitons) is 
much larger than the size of the excitons. At higher temperatures (T = 200-ftT) we observe a 
significant number of free electrons and holes (cf. Fig. |2J right panel). Due to the temperature- 
induced dissociation of the excitonic bound states a partially ionized (non-ideal) electron-hole 
plasma is formed. 

3.2. Pair distribution functions and structure factors 

This scenario is supported by the behavior of the electron-electron, hole-hole and electron-hole 
correlations pointed out in this section. The pair distribution functions and structure factors 
shown in Fig. Elare the sum over both spin projections of the particles. 

Let us first discuss the physics at short distances in terms of the pair distribution functions. 
The function g a b is the distribution of b-particles, on the average, about any a-particle. As an 
effect of the Coulomb and Fermi (statistics) repulsion, g ee and ghh are strongly depleted at small 
r. The decay of the hole-hole correlations is much stronger because of the larger masses, i.e. 
holes behave more "classically" . For the electron subsystem quantum exchange and tunnelling 
are more important and compete with the Coulomb repulsion. Furthermore the pronounced 
peak of the electron-hole pair distribution unambiguously signals the formation exciton bound 
states. This is confirmed by considering the product r 2 g e h(r) which has the physical meaning 
of the probability to find the electron at distance r away from the hole. At low temperatures, 
r2 9eh( r ) is strongly peaked around as- Increasing the temperature the excitonic peak weakens 
and finally it vanishes at about 300 K, which reflects thermal dissociation of excitons. Note that 
g ee and g^h (almost) coincide even at short distances for 200 K, i.e. the now electrons behave 
more like "classical" particles as well. Clearly all functions g a b go to unity at large distances for 
all T, which is a result of the overall uniform distributions of the charges. 
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Figure 3. (Color online) Pair distribution functions (upper panels) and structure factors (lower 
panels) at T = 50K (left panels) and T = 200K (right panels), respectively. Results are given 
for r s = 10. 



To discuss the physical behavior at large distances in more detail, it is convenient to consider 
the (charge) structure factor S a b in momentum space. According to Eq. © positive (negative) 
values of S a b indicate attraction (repulsion). From g e h the typical length scale of the attractive 
electron- hole correlations was found to be of the order of 5ae at 50 K (cf., Fig. 01 upper left 
panel). Exciton formation sets a new length scale for the structure factor as well (cf. the 
maximum in the structure factor at about 0.2 [1/ob]). In accordance with the above discussion 
this excitonic maximum is shifted to smaller k values and finally washed out if the temperature 
is raised. In the case that the maximum (minimum) of S a b is located at k = and S a b(k) 
is a monotonically decreasing (increasing) and rather structure- less function of k; the system 
closely resembles a homogeneous electron-hole plasma. Nontrivial quantum (screening) effects, 
however, come into play at low temperatures as can be seen from the wiggly behavior especially 
of the electron-electron and hole-hole charge correlations. 

3.3. Exciton density 

From the analysis of the snapshots in Sec. 3.1 we found that our complex electron-hole plasma 
contains at the same time free electrons, holes and excitons. That means, depending on the 
temperatures and particle densities, electrons and holes may appear in bound and unbound 
quantum and classical states. In a strict sense, we cannot calculate the fraction of bound and 
free states in our DPIMC simulations because of the possible overlapping of the electron "shells" , 
in particular at intermediate temperatures. Nevertheless a rough estimation of the fraction of 
the number of electron-hole bound states can be obtained by the following physical reasonings. 
The contribution of all bound and scattering states at temperature 1/(3 is given by the two- 



particle Slater sum 



oo 

S eh (r,/3) = 8vr 3 / 2 A^ £ |* a (r)| 2 exp(-/J£ a ) , (5) 

where £/ a and ^(r) are the energy (without the center of mass energy) and the wave function 
of the exciton, respectively. S e /j is, in essence, the diagonal part of the density matrix and the 
product r 2 S e / l represents the (quantum mechanical) probability for an electron and a hole to be 
separated by distance r. 

The sum over all possible states contains contributions from the discrete (d) and continuum 
(c) part of the spectrum, S e /j = S rf + S c . S rf contains all populated bound states between Eq and 
energy E' that separates bound and free (scattering) states. For low densities it is reasonable 
to choose E' ~ — ^ since higher lying bound states will be thermally ionised 12 . 

At temperatures smaller than the electron-hole binding energy and distances smaller or of 
the order of several Bohr radii, the main contribution to the Slater sum comes from electron- 
hole bound states, while at larger distances free electron-hole scattering states turn out to be 
the crucial factor. We found that the product r 2 S e /j is sharply peaked around as, just like 
r2 9eh( r ) m Fig. EJ where the maximum is most pronounced at low temperatures. Increasing the 
temperature, r 2 g e h{r) approaches the parabola r 2 related to the case of full absence of bound 
states (ideal plasma g e h{ r ) = 1)- An attractive interaction results in an increase of the pair 
distribution function at short e-h distances: g e h( r ) > 1 for < r < r b . Since the total number 
of particles is conserved this increase involves a decrease at large distances (see Fig. OJ). So 
the interval [0,r b ] with r 2 * (g e h( r ) — 1) > gives the contribution of the bound states to the 
total, bound and scattering state probability. In a strict sense, correlated e-h scattering states 
also contribute to the maximum of g e h(r), but their influence is negligible for the temperatures 
considered, where r b is of the order of several Bohr radii. For example, for T = 50-fT only 
excitons with the principle quantum number n = 1 and n = 2 are stable, which corresponds to 
a predominant electron-hole separation in the range of about 4a#. This is well reproduced by 
the procedure outlined above which yields r b ~ 3.7ob, cf. Fig. 01 (upper left part). Hence the 
fraction of the electron-hole bound states can be estimated from the ratio 
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Figure 4. (Color online) Fraction of electron-hole bound states versus temperature for different 
hole masses and dielectric constants. 



where the denominator is the normalization constant and has the physical meaning of electrons 
to be in bound or scattering states. 

Figure 0] shows the temperature dependence of the fraction of electron-hole bound states 
calculated using Eq. @ for two dielectric constants and different electron to hole mass ratios. 
Obviously, larger e's result in lower electron-hole binding energies and therefore shift the 
ionization temperature to smaller values. The same tendency is observed if m,h/m e is lowered. 
For example, from Fig. |1] we found that the ratio of the temperatures, where the fraction of 
electron- hole bound states is 50%, is the same as the corresponding ratio of the reduced electron- 
hole masses \x = m e mh/(m e + m^) (keeping e fixed). We have a ratio of about 2 for the data 
belonging to rrih = 80 and rrih = 2.1, and 1.5 for those belonging to rrih = 80 and rrih = 4.2. 
This is exactly what we would expect from the ideal Saha equation because the exciton binding 
energy is known to be proportional to the reduced electron-hole mass. 

4. Conclusions 

To summarize, we have performed a largely unbiased direct path integral quantum Monte Carlo 
simulation of a neutral non-ideal two-component Coulomb system consisting of light electrons 
and heavy holes. The careful analysis of typical many-particle configurations, pair distribution 
and charge structure functions reveals a rather comprehensive physical picture of how the system 
behaves as a function of temperature in the low-density regime. We obtained clear evidence for 
exciton formation at low temperatures and a temperature-induced dissociation of these electron- 
hole bound states. The ionization temperature strongly depends on the dielectric constant and 
the mass ratio of electrons and holes, where small values of e and large mass asymmetry stabilize 
the bound states. As found out recently, these excitonic bound states are the limiting factor for 
a possible Coulomb crystallization of the heavy holes at low densities (Hj. At large temperatures 
the system resembles an electron-hole plasma. 

Electron-hole systems with such a strong mass anisotropy are realized in intermediate valent 
Tm[Se x Tei_ x ] alloys 0. In these materials f-d hybridization provides us with a narrow dispersive 
f valence band and, as a consequence, with a large effective hole mass of the order of 50-100 
(bare) electron masses. TmSeo.45Teo.55 is, at ambient conditions, an indirect semiconductor 
having a gap of = 130 meV, where with E b x ~ 50 — 70 meV below the bottom of the d 
band an excitonic level has been observed. Applying pressure the gap can be tuned (and even 
closed) and the material is speculated to realize in the pressure region between 5 and 11 kbar an 
excitonic insulator, the search for which has been run for a long time. A necessary precondition 
is the existence of a large number of (up to ~ 10 20 ) excitons with intermediate size (in order 
to avoid too strong overlap of the excitonic bound states) [S]. Surprisingly, in TmSeo.45Teo.55 
the excitonic phase then is predicted to be stable at rather high temperatures (up to 200 K). 
Using the Tm[Se x Tei_ x ] parameters within our (surely oversimplified) two-component Coulomb 
model DPIMC simulation, we could at least corroborate this belief. We found that (i) the 
fraction of excitonic bound states amounts to 80-90 %, (ii) the bound states will be stable up to 
100-150 K, and (iii) the mass asymmetry between holes and electron is crucial therefore. As yet 
we are not in the position to detect the condensed excitonic phase. Feasible, however, should be 
the investigation of the density dependence of the various quantities. Here we expect to see the 
formation of bi-exciton bound states and finally the Mott transition at large particle densities n a 
(i.e. small r s ). Note that the excitonic-insulator semi-metal (Mott) transition has been detected 
to occur in TmSeo.45Teo.55 at a pressure of about 13-14 kbar for T — > 0. Work along this line is 
in progress. 
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